# Replication code for: 
    # Are Democratic Innovations Legitimate to the Broader Public? Evidence from Survey Experiments in Ghana
    # Eric Kramon
    # British Journal of Political Science

# create "tables" and "figs" in the current working directory
dir.create("tables", showWarnings = FALSE)
dir.create("figs",    showWarnings = FALSE)


# ----  Packages ----
library(tidyverse)
library(ggplot2)
library(broom)
library(estimatr)
library(patchwork)
library(forcats)
library(scales)
library(lemon)
library(xtable)
library(fastDummies)
library(knitr)
library(kableExtra)
library(nnet)        # multinom()
library(texreg)      # regression tables


# ---- Bring in Data ----

# Individual level data (used for descriptive statistics)
dat <- readRDS("dat_rep.rds")

# Conjoint 1 data -- row is the participant-profile dyad
dc <- readRDS("conjoint1_rep.rds")

# Binary outcomes for Democratic and Fair variables (created once here)
dc <- dc %>% 
  mutate(
    democratic_binary = ifelse(democratic %in% c("Extremely democratic","Very democratic","Democratic"), 1, 0),
    fair_binary       = ifelse(fair %in% c("Extremely fair","Very fair","Fair"), 1, 0)
  )


# Conjoint 2 data
c2L <- readRDS("conjoint2_rep.rds")


# Labels used across plots
labs  <- c("Status quo", "Citizen-Elite", "Participatory", "Mini-Public")
labs2 <- c("Citizen-Elite", "Participatory", "Mini-Public")


# =============== Main Text Figures ===============

# ---------- Figure 1 (Satisfaction with Democracy in Ghana over Time) ----------
## values inputted by hand from Afrobarometer website

response <- c(
  rep("Not at all", 10),
  rep("Not very", 10),
  rep("Fairly", 10),
  rep("Very", 10)
)
prop <- c(
  15.6, 6.3, 4.3, 6.4, 5.7, 15.3, 5.1, 10.4, 17.4, 21.0,
  16.4, 10.1, 10.9, 9.1, 19.3, 20.1, 15.5, 21.1, 31.3, 28.3,
  36.6, 26.5, 32.4, 25.6, 48.7, 35.4, 45.5, 40.5, 31.1, 30.4,
  17.3, 19.6, 37.7, 54.4, 25.1, 20.1, 32.6, 25.8, 19.4, 19.5
)
year <- rep(c(1999, 2002, 2005, 2008, 2012, 2014, 2017, 2019, 2022, 2024), 4)

s <- tibble(response, prop, year) %>%
  mutate(satisfied = ifelse(response %in% c("Not at all","Not very"), "Not Satisfied", "Satisfied"))

plotData <- s %>%
  group_by(year, satisfied) %>%
  summarize(prop = sum(prop)/100, .groups = "drop") %>%
  mutate(var = prop*(1-prop), se = sqrt(var/1000),
         upper = prop + 1.96*se, lower = prop - 1.96*se)

p_fig_satisfaction <- ggplot(plotData, aes(x = year, y = prop, color = satisfied, shape = satisfied, ymax = upper, ymin = lower)) +
  geom_line(linetype = 2) +
  geom_point() +
  geom_errorbar(width = .1) +
  theme_bw() +
  scale_color_viridis_d(option = "inferno", begin = .2, end = .7) +
  labs(y = "Proportion", x = "Year", color = "", shape = "", caption = "Data Source: Afrobarometer Surveys") +
  theme(legend.position="bottom")

p_fig_satisfaction

ggsave("figs/figure1.pdf", p_fig_satisfaction, width = 6, height = 4)



# ---------- Figure: 2 (Full Conjoint AMCE on Binary Choice) ----------
dc <- dc %>% 
  mutate(pr2 = fct_relevel(as.factor(process), "1"), 
         out2 = fct_relevel(as.factor(outcome), "2"),
         endorse = as.factor(endorse))

res <- lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = dc) %>% 
  tidy() %>% 
  mutate(Attribute = case_when(
           term %in% c("pr22","pr23","pr24") ~ "Process",
           term %in% c("endorse2","endorse3") ~ "Endorsement",
           term == "out21" ~ "Outcome",
           TRUE ~ NA_character_
         )) %>% 
  select(term, estimate, conf.low, conf.high, Attribute) %>% 
  filter(term != "(Intercept)")

p1 <- tibble(term = c("Status Quo","No Endorsement","No Project"),
             estimate = 0, conf.low = 0, conf.high = 0,
             Attribute = c("Process","Endorsement","Outcome"))
pd <- bind_rows(res, p1)

pd1 <- pd %>% filter(Attribute == "Process")
pd2 <- pd %>% filter(Attribute == "Endorsement")
pd3 <- pd %>% filter(Attribute == "Outcome")

f1 <- ggplot(pd1, aes(y = estimate, x = term, ymin = conf.low, ymax = conf.high)) +
  geom_point() + geom_errorbar(width = 0) + theme_bw() + coord_flip() +
  ylim(-.1, .3) + geom_hline(yintercept = 0, linetype = 2, color = "darkgrey") +
  scale_x_discrete(labels = c("Citizen-Elite","Participatory","Mini-Public","Status Quo")) +
  facet_grid(rows = vars(Attribute), switch = "y") + labs(y = "", x = "") +
  theme(strip.placement = "outside", strip.text.y = element_text(angle = 0))

f2 <- ggplot(pd2, aes(y = estimate, x = term, ymin = conf.low, ymax = conf.high)) +
  geom_point() + geom_errorbar(width = 0) + theme_bw() + coord_flip() +
  ylim(-.1, .3) + geom_hline(yintercept = 0, linetype = 2, color = "darkgrey") +
  scale_x_discrete(labels = c("NPP Endorsement","NDC Endorsement","No Endorsement")) +
  facet_grid(rows = vars(Attribute), switch = "y") + labs(y = "", x = "") +
  theme(strip.placement = "outside", strip.text.y = element_text(angle = 0))

f3 <- ggplot(pd3, aes(y = estimate, x = term, ymin = conf.low, ymax = conf.high)) +
  geom_point() + geom_errorbar(width = 0) + theme_bw() + coord_flip() +
  ylim(-.1, .3) + geom_hline(yintercept = 0, linetype = 2, color = "darkgrey") +
  scale_x_discrete(labels = c("No Project","Project")) +
  facet_grid(rows = vars(Attribute), switch = "y") +
  labs(y = "Average Marginal Component Effect", x = "") +
  theme(strip.placement = "outside", strip.text.y = element_text(angle = 0))

p_fig_main <- f1 / f2 / f3
p_fig_main

ggsave("figs/figure2.pdf", p_fig_main, width = 6, height = 8)


# ---------- Figure 3  (Legitimacy Evaluations) ----------

pres2 <- dc %>% group_by(process) %>% 
  do(tidy(lm_robust(as.numeric(fair) ~ 1, data = ., cluster = ID))) %>% 
  mutate(var = "Fairness")

pres3 <- dc %>% 
  group_by(process) %>% 
  do(tidy(lm_robust(as.numeric(democratic) ~ 1, data = ., cluster = ID))) %>% 
  mutate(var = "Democratic")

e1 <- bind_rows(pres2, pres3)

p_fig_process <- ggplot(e1, aes(y=estimate, x=factor(process),
                                ymin = conf.low, ymax = conf.high,
                                color = var, shape = var)) +
  geom_point(position=position_dodge(width = -.6)) +
  geom_errorbar(width = 0, position=position_dodge(width = -.6)) +
  theme_bw() +
  scale_x_discrete(breaks=seq(1,4,1), labels = labs) +
  labs(x = "Process", y = "", color = "", shape = "")  +
  scale_y_symmetric(mid = 4) +
  geom_hline(yintercept = 4, linetype = 2, color = "grey") +
  scale_colour_manual(values = c("#C67800", "#205C8A")) +
  theme(legend.position="bottom")

p_fig_process
ggsave("figs/figure3.pdf", p_fig_process, width = 6, height = 4)

# ---------- Figure: 4 (Preferences by Outcome) ----------
inter <- dc %>% mutate(project = as.factor(outcome)) %>%
  group_by(process, project) %>% do(tidy(lm_robust(choice ~ 1, data = ., cluster = ID)))

p_fig_outcome1 <- ggplot(inter, aes(y=estimate, x=process, color=project, shape=project)) +
  geom_point(position=position_dodge(width = -.4)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position=position_dodge(width = -.4)) +
  theme_bw() + labs(y = "Probability of Selection", x = "Process", color = "", shape = "") +
  scale_x_discrete(breaks=seq(1,4,1), labels = labs) +
  scale_colour_manual(name=NULL, labels=c("Project","No Project"), values = c("#C67800","#205C8A")) +
  scale_shape_manual(name=NULL, labels=c("Project","No Project"), values = c(19,17)) +
  scale_y_symmetric(mid = .5) + geom_hline(yintercept = .5, linetype = 2, color = "grey") +
  theme(legend.position="bottom")


p_fig_outcome1

ggsave("figs/figure4.pdf", p_fig_outcome1, width = 6, height = 4.5)


# ---------- Figure 5 (Legitimacy by Outcome) ----------

inter_fair <- dc %>% mutate(project = as.factor(outcome)) %>%
  group_by(project, process) %>% do(tidy(lm_robust(as.numeric(fair) ~ 1, data = ., cluster = ID)))

fa <- ggplot(inter_fair, aes(y=estimate, x=factor(process), color=project, shape=project)) +
  geom_point(position=position_dodge(width = -.6)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position=position_dodge(width = -.6)) +
  theme_bw() + xlab(" ") + ylab(" ") + scale_x_discrete(breaks=seq(1,4,1), labels=labs) +
  ggtitle("(a) Fairness") +
  scale_colour_manual(values = c("#C67800", "#205C8A")) +
  scale_shape_manual(values = c(19,17)) +
  scale_y_symmetric(mid = 4) + geom_hline(yintercept = 4, linetype = 2, color = "grey") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

inter_dem <- dc %>% mutate(project = as.factor(outcome)) %>%
  group_by(project, process) %>% do(tidy(lm_robust(as.numeric(democratic) ~ 1, data = ., cluster = ID)))
fb <- ggplot(inter_dem, aes(y=estimate, x=factor(process), color=project, shape=project)) +
  geom_point(position=position_dodge(width = -.6)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position=position_dodge(width = -.6)) +
  theme_bw() + xlab(" ") + ylab(" ") + scale_x_discrete(breaks=seq(1,4,1), labels=labs) +
  ggtitle("(b) Democratic") +
  scale_colour_manual(values = c("#C67800", "#205C8A")) +
  scale_shape_manual(values = c(19,17)) +
  scale_y_symmetric(mid = 4) + geom_hline(yintercept = 4, linetype = 2, color = "grey") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p_fig_outcome2 <- fa + fb + plot_layout(guides="collect", ncol = 2) & theme(legend.position = 'bottom')
p_fig_outcome2

ggsave("figs/figure5.pdf", p_fig_outcome2, width = 7, height = 4.5)

# ---------- Figure 6(Experiment 2 — Preferences on Reforms) ----------

c2L <- c2L %>% mutate(effective01 = rescale(as.numeric(effective)))

choice <- c2L %>% group_by(reform) %>% 
  do(tidy(lm_robust(choice ~ 1, data = ., cluster = ID))) %>% 
  mutate(var = "More Effective (choice)")

dem    <- c2L %>% group_by(reform) %>% 
  do(tidy(lm_robust(democratic ~ 1, data = ., cluster = ID))) %>% 
  mutate(var = "More Democratic (choice)")

dp <- bind_rows(choice, dem)

p_fig_main2 <- ggplot(dp, aes(y=estimate, x=reform, color = var, shape = var)) +
  geom_point(position=position_dodge(width = -.6)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position=position_dodge(width = -.6)) +
  theme_bw() + labs(color = "", y = "Probability of Selection", x = "", shape = "") +
  scale_x_discrete(breaks=seq(1,4,1), labels=c("Elections","Direct Democracy","Citizen-Elite","Mini-Publics")) +
  scale_colour_manual(values = c("#C67800","#205C8A")) +
  theme(legend.position="bottom") +
  scale_y_symmetric(mid = .5) + geom_hline(yintercept = .5, linetype = 2, color = "grey")

p_fig_main2
ggsave("figs/figure6.pdf", p_fig_main2, width = 6, height = 4.5)

# =============== Main Text Table ===============

# ---------- Table 1 (Descriptive Statistics) ----------
dat1 <- fastDummies::dummy_cols(dat, select_columns = c("age1","ed1"))
n <- nrow(dat1)
summ <- dat1 %>% summarize(
  Female = mean(female, na.rm = TRUE),
  Rural = mean(rural, na.rm = TRUE),
  `Age 18-24` = mean(`age1_16-24 years old`, na.rm = TRUE),
  `Age 25-34` = mean(`age1_25-34 years old`, na.rm = TRUE),
  `Age 35-44` = mean(`age1_35-44 years old`, na.rm = TRUE),
  `Age 45-54` = mean(`age1_45-54 years old`, na.rm = TRUE),
  `Age 55-64` = mean(`age1_55-64 years old`, na.rm = TRUE),
  `Age 65+`   = mean(`age1_65+ years old`, na.rm = TRUE),
  `No Formal Education` = mean(ed1_0, na.rm = TRUE),
  `Primary Education`   = mean(ed1_1, na.rm = TRUE),
  `Secondary Education` = mean(ed1_2, na.rm = TRUE),
  `University and above Education` = mean(ed1_3, na.rm = TRUE),
  `NPP Supporter` = sum(npp_partisan, na.rm = TRUE)/n,
  `NDC Supporter` = sum(ndc_partisan, na.rm = TRUE)/n
)
summ2 <- tidyr::pivot_longer(summ, cols = everything()) %>% rename(`Survey Mean` = value, Variable = name)
ab <- c(0.50, 0.41, 0.23, 0.31, 0.21, 0.12, NA , NA , 0.15, 0.25, 0.41, 0.19, 0.23, 0.19)
summ3 <- bind_cols(summ2, tibble(`Afrobarometer Mean` = ab))

tab1 <- xtable(summ3, digits = 3)
tab1
print(tab1, file = "tables/tbl-descriptives.tex", comment = FALSE, floating = FALSE, include.rownames = FALSE)

# =============== Online Supplementary Material ===============

# ---------- Figure C1 ----------
dc <- dc %>% mutate(pr2 = fct_relevel(process, "1"), out2 = fct_relevel(outcome, "2"))
res1 <- dc %>% filter(round == 1) %>% lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = .) %>% tidy() %>%
  mutate(Attribute = case_when(term %in% c("pr22","pr23","pr24") ~ "Process", term %in% c("endorse2","endorse3") ~ "Endorsement", term == "out21" ~ "Outcome")) %>%
  select(term, estimate, conf.low, conf.high, Attribute) %>% filter(term != "(Intercept)")
res2 <- dc %>% filter(round == 2) %>% lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = .) %>% tidy() %>%
  mutate(Attribute = case_when(term %in% c("pr22","pr23","pr24") ~ "Process", term %in% c("endorse2","endorse3") ~ "Endorsement", term == "out21" ~ "Outcome")) %>%
  select(term, estimate, conf.low, conf.high, Attribute) %>% filter(term != "(Intercept)")
p1 <- tibble(term = c("Status Quo","No Endorsement","No Project"), estimate = 0, conf.low = 0, conf.high = 0, Attribute = c("Process","Endorsement","Outcome"))
pdr1 <- bind_rows(res1, p1) %>% mutate(Round = 1)
pdr2 <- bind_rows(res2, p1) %>% mutate(Round = 2)
pd <- bind_rows(pdr1, pdr2)

pd1 <- pd %>% filter(Attribute == "Process")
pd2 <- pd %>% filter(Attribute == "Endorsement")
pd3 <- pd %>% filter(Attribute == "Outcome")

g1 <- ggplot(pd1, aes(y=estimate, x=term, ymin=conf.low, ymax=conf.high, color=factor(Round), shape=factor(Round))) +
  geom_point(position=position_dodge(.25)) + geom_errorbar(position=position_dodge(.25), width=0) +
  theme_bw() + coord_flip() + ylim(-.1, .35) + geom_hline(yintercept=0, linetype=2, color="darkgrey") +
  scale_x_discrete(labels=c("Citizen-Elite","Participatory","Mini-Public","Status Quo")) +
  facet_grid(rows = vars(Attribute), switch="y") + labs(y="", x="", shape="Round", color="Round") +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom")

g2 <- ggplot(pd2, aes(y=estimate, x=term, ymin=conf.low, ymax=conf.high, color=factor(Round), shape=factor(Round))) +
  geom_point(position=position_dodge(.25)) + geom_errorbar(position=position_dodge(.25), width=0) +
  theme_bw() + coord_flip() + ylim(-.1, .35) + geom_hline(yintercept=0, linetype=2, color="darkgrey") +
  scale_x_discrete(labels=c("NPP Endorsement","NDC Endorsement","No Endorsement")) +
  facet_grid(rows = vars(Attribute), switch="y") + labs(y="", x="", color="Round", shape="Round") +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom")

g3 <- ggplot(pd3, aes(y=estimate, x=term, ymin=conf.low, ymax=conf.high, color=factor(Round), shape=factor(Round))) +
  geom_point(position=position_dodge(.25)) + geom_errorbar(position=position_dodge(.25), width=0) +
  theme_bw() + coord_flip() + ylim(-.1, .35) + geom_hline(yintercept=0, linetype=2, color="darkgrey") +
  scale_x_discrete(labels=c("No Project","Project")) +
  facet_grid(rows = vars(Attribute), switch="y") + labs(y="Average Marginal Component Effect", x="", color="Round", shape="Round") +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom")

p_fig_mainOrder <- g1 / g2 / g3 + plot_layout(guides="collect") & theme(legend.position="bottom")

p_fig_mainOrder
ggsave("figs/fig-mainOrder.pdf", p_fig_mainOrder, width = 6, height = 8)

# ---------- Figure C2 (E1 profile-order test) ----------
res1 <- dc %>% filter(order == 1) %>% lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = .) %>% tidy() %>%
  mutate(Attribute = case_when(term %in% c("pr22","pr23","pr24") ~ "Process", term %in% c("endorse2","endorse3") ~ "Endorsement", term == "out21" ~ "Outcome")) %>%
  select(term, estimate, conf.low, conf.high, Attribute) %>% filter(term != "(Intercept)")
res2 <- dc %>% filter(order == 2) %>% lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = .) %>% tidy() %>%
  mutate(Attribute = case_when(term %in% c("pr22","pr23","pr24") ~ "Process", term %in% c("endorse2","endorse3") ~ "Endorsement", term == "out21" ~ "Outcome")) %>%
  select(term, estimate, conf.low, conf.high, Attribute) %>% filter(term != "(Intercept)")
p1 <- tibble(term = c("Status Quo","No Endorsement","No Project"), estimate = 0, conf.low = 0, conf.high = 0, Attribute = c("Process","Endorsement","Outcome"))
pdr1 <- bind_rows(res1, p1) %>% mutate(Round = 1)
pdr2 <- bind_rows(res2, p1) %>% mutate(Round = 2)
pd <- bind_rows(pdr1, pdr2)
pd1 <- pd %>% filter(Attribute == "Process")
pd2 <- pd %>% filter(Attribute == "Endorsement")
pd3 <- pd %>% filter(Attribute == "Outcome")

g1 <- ggplot(pd1, aes(y=estimate, x=term, ymin=conf.low, ymax=conf.high, color=factor(Round), shape=factor(Round))) +
  geom_point(position=position_dodge(.25)) + geom_errorbar(position=position_dodge(.25), width=0) + theme_bw() + coord_flip() +
  ylim(-.1, .35) + geom_hline(yintercept=0, linetype=2, color="darkgrey") +
  scale_x_discrete(labels=c("Citizen-Elite","Participatory","Mini-Public","Status Quo")) +
  facet_grid(rows=vars(Attribute), switch="y") + labs(y="", x="", shape="Order", color="Order") +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom")

g2 <- ggplot(pd2, aes(y=estimate, x=term, ymin=conf.low, ymax=conf.high, color=factor(Round), shape=factor(Round))) +
  geom_point(position=position_dodge(.25)) + geom_errorbar(position=position_dodge(.25), width=0) + theme_bw() + coord_flip() +
  ylim(-.1, .35) + geom_hline(yintercept=0, linetype=2, color="darkgrey") +
  scale_x_discrete(labels=c("NPP Endorsement","NDC Endorsement","No Endorsement")) +
  facet_grid(rows=vars(Attribute), switch="y") + labs(y="", x="", color="Order", shape="Order") +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom")

g3 <- ggplot(pd3, aes(y=estimate, x=term, ymin=conf.low, ymax=conf.high, color=factor(Round), shape=factor(Round))) +
  geom_point(position=position_dodge(.25)) + geom_errorbar(position=position_dodge(.25), width=0) + theme_bw() + coord_flip() +
  ylim(-.1, .35) + geom_hline(yintercept=0, linetype=2, color="darkgrey") +
  scale_x_discrete(labels=c("No Project","Project")) +
  facet_grid(rows=vars(Attribute), switch="y") + labs(y="Average Marginal Component Effect", x="", color="Order", shape="Order") +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom")

p_fig_mainProfileOrder <- g1 / g2 / g3 + plot_layout(guides="collect") & theme(legend.position="bottom")

p_fig_mainProfileOrder

ggsave("figs/fig-mainProfileOrder.pdf", p_fig_mainProfileOrder, width = 6, height = 8)


# ---------- Figure C3 (Experiment 2 profile-order check) ----------
c2L <- c2L %>% mutate(effective01 = rescale(as.numeric(effective)))

# first profile listed
choice1 <- c2L %>% filter(order == 1) %>% group_by(reform) %>% do(tidy(lm_robust(choice ~ 1, data = ., cluster = ID))) %>% mutate(var = "More Effective (choice)")
dem1    <- c2L %>% filter(order == 1) %>% group_by(reform) %>% do(tidy(lm_robust(democratic ~ 1, data = ., cluster = ID))) %>% mutate(var = "More Democratic (choice)")
dp1 <- bind_rows(choice1, dem1)
p1 <- ggplot(dp1, aes(y=estimate, x=reform, color = var, shape = var)) +
  geom_point(position=position_dodge(width = -.6)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position=position_dodge(width = -.6)) +
  theme_bw() + labs(title="(a) First Profile Listed", color="", y="Probability of Selection", x="", shape="") +
  scale_x_discrete(breaks=seq(1,4,1), labels=c("Elections","Direct Democracy","Citizen-Elite","Mini-Publics")) +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom") +
  scale_y_symmetric(mid = .5) + geom_hline(yintercept=.5, linetype=2, color="grey")

# second profile listed
choice2 <- c2L %>% filter(order == 2) %>% group_by(reform) %>% do(tidy(lm_robust(choice ~ 1, data = ., cluster = ID))) %>% mutate(var = "More Effective (choice)")
dem2    <- c2L %>% filter(order == 2) %>% group_by(reform) %>% do(tidy(lm_robust(democratic ~ 1, data = ., cluster = ID))) %>% mutate(var = "More Democratic (choice)")
dp2 <- bind_rows(choice2, dem2)
p2 <- ggplot(dp2, aes(y=estimate, x=reform, color = var, shape = var)) +
  geom_point(position=position_dodge(width = -.6)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position=position_dodge(width = -.6)) +
  theme_bw() + labs(title="(b) Second Profile Listed", color="", y="Probability of Selection", x="", shape="",
                    caption = "Dependent variables are the two binary measures. 95% CIs clustered on subject.") +
  scale_x_discrete(breaks=seq(1,4,1), labels=c("Elections","Direct Democracy","Citizen-Elite","Mini-Publics")) +
  scale_colour_manual(values=c("#C67800","#205C8A")) + theme(legend.position="bottom") +
  scale_y_symmetric(mid = .5) + geom_hline(yintercept=.5, linetype=2, color="grey")

p_fig_exp2carry <- p1 / p2

p_fig_exp2carry
ggsave("figs/fig-exp2carry.pdf", p_fig_exp2carry, width = 6.5, height = 8)



# ---------- Figure E1 (Experiment 2 effectiveness 1–7 scale) ----------
eff <- c2L %>% group_by(reform) %>% do(tidy(lm_robust(as.numeric(effective) ~ 1, data = ., cluster = ID))) %>% mutate(var = "Effectiveness (scale)")
p_fig_main2eff <- ggplot(eff, aes(y=estimate, x= fct_reorder(reform, -estimate))) +
  geom_point(color = "steelblue4") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, color = "steelblue4") +
  theme_bw() + labs(y="How Effective (1–7)", x="") +
  scale_x_discrete(breaks=seq(1,4,1), labels=c("Elections","Direct Democracy","Citizen-Elite","Mini-Publics"))

p_fig_main2eff
ggsave("figs/fig-main2eff.pdf", p_fig_main2eff, width = 6, height = 4.5)

# =============== Appendix Regression Tables  ===============

# ---------- Table D1 Full regression results from Exp 1 ----------
res_all <- lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = dc)
res_ndc <- dc %>% filter(ndc_partisan == 1) %>% lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = .)
res_npp <- dc %>% filter(npp_partisan == 1) %>% lm_robust(choice ~ pr2 + endorse + out2, cluster = ID, data = .)
texreg(list(res_all, res_ndc, res_npp),
       stars = c(0.10, 0.05, 0.01),
       digits = 3,
       include.ci = FALSE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE,
       include.f = FALSE,
       custom.model.names = c("All Respondents", "NDC Supporters", "NPP Supporters"),
       custom.coef.names = c("(Intercept)", "Citizen-Elite", "Participatory", "Mini-Public", "NPP Endorsement", "NDC Endorsement", "Receives Project"),
       custom.gof.names = c("N (respondent-process)", "N (respondents)"),
       float.pos = "",
       custom.note = "\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-regressionAll.tex"
)

# ---------- Table D2 (E1: Citizen-Elite as reference; choice, democratic, fair) ----------
dc <- dc %>% mutate(
  pr_f = case_match(as.numeric(pr2), 1 ~ "Status Quo", 2 ~ "Citizen-Elite", 3 ~ "Participatory", 4 ~ "Mini-Public"),
  pr_f = factor(pr_f)
)
dc$pr_f <- relevel(dc$pr_f, ref = "Citizen-Elite")
res_choice <- lm_robust(choice ~ pr_f + endorse + out2, cluster = ID, data = dc)
res_dem    <- lm_robust(as.numeric(democratic) ~ pr_f + endorse + out2, cluster = ID, data = dc)
res_fair   <- lm_robust(as.numeric(fair) ~ pr_f + endorse + out2, cluster = ID, data = dc)
texreg(list(res_choice, res_dem, res_fair),
       stars = c(0.10, 0.05, 0.01),
       digits = 3,
       include.ci = FALSE,
       include.rsquared = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE,
       include.f = FALSE,
       custom.model.names = c("Choice", "Democratic", "Fairness"),
       custom.coef.names = c("(Intercept)", "Mini-Public", "Participatory", "Status Quo", "NPP Endorsement", "NDC Endorsement", "Receives Project"),
       custom.gof.names = c("N (respondent-process)", "N (respondents)"),
       float.pos = "",
       custom.note = "\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-regressionCEbaseline.tex"
)

# ---------- Table B1 (E1 balance: predictors of assignment to process) ----------
dc_bal <- dc %>% mutate(
  prf = as.numeric(pr2),
  pr_fac = case_match(prf, 1 ~ "Status Quo", 2 ~ "Citizen-Elite", 3 ~ "Participatory", 4 ~ "Mini-Public"),
  Female = female, Rural = rural, Age = age, Education = education,
  `NPP Partisan` = npp_partisan, `NDC Partisan` = ndc_partisan
)
process_mod <- nnet::multinom(pr_fac ~ Female + Rural + Age + Education + `NPP Partisan` + `NDC Partisan`, data = dc_bal, trace = FALSE)
tab_outcomeProcess <- tidy(process_mod) %>% 
  select(y.level, term, estimate, std.error, p.value) %>%
  rename(Outcome = y.level, Predictor = term, Estimate = estimate, `Standard Error` = std.error, `P-Value` = p.value) %>%
  mutate(Predictor = gsub(":", " x ", Predictor)) %>%
  knitr::kable(format = "latex", booktabs = TRUE, digits = 3)
writeLines(tab_outcomeProcess, "tables/tbl-outcomeProcess.tex")

# ---------- Table B2 (E1 balance: predictors of endorsement) ----------
dc_end <- dc %>% mutate(
  endorsef = as.numeric(endorse),
  endorse_fac = case_match(endorsef, 1 ~ "No Endorsement", 2 ~ "NPP Endorsement", 3 ~ "NDC Endorsement"),
  Female = female, Rural = rural, Age = age, Education = education,
  `NPP Partisan` = npp_partisan, `NDC Partisan` = ndc_partisan
)
endorse_mod <- nnet::multinom(endorse_fac ~ Female + Rural + Age + Education + `NPP Partisan` + `NDC Partisan`, data = dc_end, trace = FALSE)
tab_endorseBalance <- tidy(endorse_mod) %>% 
  select(y.level, term, estimate, std.error, p.value) %>%
  rename(Outcome = y.level, Predictor = term, Estimate = estimate, `Standard Error` = std.error, `P-Value` = p.value) %>%
  mutate(Predictor = gsub(":", " x ", Predictor)) %>%
  knitr::kable(format = "latex", booktabs = TRUE, digits = 3)
writeLines(tab_endorseBalance, "tables/tbl-endorseBalance.tex")

# ---------- Table B3 (E1 balance: predictors of receiving project) ----------
dc_out <- dc %>% 
  mutate(o = ifelse(out2 == "1", 1, 0),
         Female = female, Rural = rural, Age = age, Education = education,
         `NPP Partisan` = npp_partisan, `NDC Partisan` = ndc_partisan) 
out_mod <- lm(o ~ Female + Rural + Age + Education + `NPP Partisan` + `NDC Partisan`, data = dc_out)
texreg(list(out_mod),
       stars = c(0.10, 0.05, 0.01),
       digits = 3, include.ci = FALSE, include.rsquared = FALSE, include.adjrs = FALSE,
       include.rmse = FALSE, include.f = FALSE, float.pos = "",
       custom.note = "\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Standard errors in parentheses.}",
       file = "tables/tbl-outcomeBalance.tex")

# ---------- Table B4 (E2 balance: predictors of assignment to reform) ----------
c2L_bal <- c2L %>% mutate(
  ref2 = as.numeric(reform),
  pr_fac = case_match(ref2, 1 ~ "Direct Democracy", 2 ~ "Citizen-Elite", 3 ~ "Elections", 4 ~ "Mini-Publics"),
  Female = female, Rural = rural, Age = age, Education = education,
  `NPP Partisan` = npp_partisan, `NDC Partisan` = ndc_partisan
)
process2_mod <- nnet::multinom(pr_fac ~ Female + Rural + Age + Education + `NPP Partisan` + `NDC Partisan`, data = c2L_bal, trace = FALSE)
tab_exp2balance <- tidy(process2_mod) %>% 
  select(y.level, term, estimate, std.error, p.value) %>%
  rename(Outcome = y.level, Predictor = term, Estimate = estimate, `Standard Error` = std.error, `P-Value` = p.value) %>%
  mutate(Predictor = gsub(":", " x ", Predictor)) %>%
  knitr::kable(format = "latex", booktabs = TRUE, digits = 3)
writeLines(tab_exp2balance, "tables/tbl-exp2balance.tex")

# ---------- Table D3 (E1) ----------
m_rural <- lm_robust(choice ~ pr2*rural + endorse + out2, cluster = ID, data = dc)
texreg(list(m_rural), stars = c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names = c(""), custom.coef.names = c("(Intercept)", "Citizen-Elite", "Participatory", "Mini-Public", "Rural", "NPP Endorsement", "NDC Endorsement", "Receives Project", "Rural x Citizen-Elite", "Rural x Participatory", "Rural x Mini-Public"),
       custom.gof.names = c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-rural.tex")

# ---------- Table D4 (E1) ----------

m_age <- lm_robust(choice ~ pr2*age + endorse + out2, cluster = ID, data = dc)
texreg(list(m_age), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Citizen-Elite", "Participatory", "Mini-Public", "Age", "NPP Endorsement", "NDC Endorsement", "Receives Project", "Age x Citizen-Elite", "Age x Participatory", "Age x Mini-Public"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-age.tex")

# ---------- Table 5(E1) ----------

m_edu <- lm_robust(choice ~ pr2*education + endorse + out2, cluster = ID, data = dc)
texreg(list(m_edu), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Citizen-Elite", "Participatory", "Mini-Public", "Education", "NPP Endorsement", "NDC Endorsement", "Receives Project", "Education x Citizen-Elite", "Education x Participatory", "Education x Mini-Public"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-education.tex")


# ---------- Table D6 (E1) ----------

m_npp <- lm_robust(choice ~ pr2*npp_partisan + endorse + out2, cluster = ID, data = dc)
texreg(list(m_npp), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Citizen-Elite", "Participatory", "Mini-Public", "NPP Partisan", "NPP Endorsement", "NDC Endorsement", "Receives Project", "NPP Partisan x Citizen-Elite", "NPP Partisan x Participatory", "NPP Partisan x Mini-Public"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-npp.tex")

# ---------- Table D7 (E1) ----------

m_ndc <- lm_robust(choice ~ pr2*ndc_partisan + endorse + out2, cluster = ID, data = dc)
texreg(list(m_ndc), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Citizen-Elite", "Participatory", "Mini-Public", "NDC Partisan", "NPP Endorsement", "NDC Endorsement", "Receives Project", "NDC Partisan x Citizen-Elite", "NDC Partisan x Participatory", "NDC Partisan x Mini-Public"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-ndc.tex")

# ---------- Table E1 ----------
m_rural_e2 <- lm_robust(choice ~ reform*rural, data = c2L, cluster = ID)
texreg(list(m_rural_e2), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Direct Democracy", "Citizen-Elite", "Mini-Publics", "Rural", "Rural x Direct Democracy", "Rural x Citizen-Elite", "Rural x Mini-Publics"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-rurale2.tex")


# ---------- Table E2 ----------

m_age_e2 <- lm_robust(choice ~ reform*age, data = c2L, cluster = ID)
texreg(list(m_age_e2), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Direct Democracy", "Citizen-Elite", "Mini-Publics", "Age", "Age x Direct Democracy", "Age x Citizen-Elite", "Age x Mini-Publics"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-agee2.tex")

# ---------- Table E3 ----------

m_edu_e2 <- lm_robust(choice ~ reform*education, data = c2L, cluster = ID)
texreg(list(m_edu_e2), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Direct Democracy", "Citizen-Elite", "Mini-Publics", "Education", "Education x Direct Democracy", "Education x Citizen-Elite", "Education x Mini-Publics"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-educatione2.tex")

# ---------- Table E4 ----------

m_npp_e2 <- lm_robust(choice ~ reform*npp_partisan, data = c2L, cluster = ID)
texreg(list(m_npp_e2), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Direct Democracy", "Citizen-Elite", "Mini-Publics", "NPP Partisan", "NPP Partisan x Direct Democracy", "NPP Partisan x Citizen-Elite", "NPP Partisan x Mini-Publics"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-nppe2.tex")

# ---------- Table E5 ----------

m_ndc_e2 <- lm_robust(choice ~ reform*ndc_partisan, data = c2L, cluster = ID)
texreg(list(m_ndc_e2), stars=c(0.10,0.05,0.01), digits=3, include.ci=FALSE,
       include.rsquared=FALSE, include.adjrs=FALSE, include.rmse=FALSE, include.f=FALSE,
       custom.model.names=c(""), custom.coef.names=c("(Intercept)", "Direct Democracy", "Citizen-Elite", "Mini-Publics", "NDC Partisan", "NDC Partisan x Direct Democracy", "NDC Partisan x Citizen-Elite", "NDC Partisan x Mini-Publics"),
       custom.gof.names=c("N (respondent-process)", "N (respondents)"), float.pos="",
       custom.note="\\parbox{.95\\linewidth}{\\vspace{6pt} \\normalsize %stars. Robust standard errors clustered by respondent in parentheses.}",
       file = "tables/tbl-ndce2.tex")

